% Based on replication codes for Dario Caldara and Christophe Kamps (2017), The Analytics of SVARs: A
% Unified Framework to Measure Fiscal Multipliers, Review of Economic Studies (2017) 84, 1015–1040
%
% Copyright: 2017 Dario Caldara and Christophe Kamps
% Copyright: 2020-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer

clear all
rng("default")
%-----------------
% Housekeeping
%----------------
addpath('./auxfiles')
if ~isfolder('results')
    mkdir('.','results');
end
if ~isfolder('Figures')
    mkdir('.','Figures');
end


options = optimoptions('fsolve','Display','none');
%--------------------
% Model Specification
%--------------------
p               = 4;      % Number of lags
nex_            = 1;      % Constant
nP              = 2; % Number of fiscal variables
T0              = p;      % Size of pre-sample for Minnesota Prior
Horizon         = 41;     % Horizon for Impulse Responses
nd              = 20000;  % Number of draws in MC chain
bburn           = 0*nd; % Share of draws to burn; 0 here due to conjugate prior not inducing a Markov Chain

ptileVEC = [0.05 0.16 0.50 0.84 0.95]; % Percentiles for posterior distributions

do_BP        = 1;   % Execute Blanchard Perotti (2002) identification
do_PF        = 1;   % Execute penalty function identification

% Number of periods for penalty function identification
nperPF = 1;


%--------------------
% Load Data
%--------------------

[data_US, text]=xlsread("US_data.xlsx");
data.var_names={'G','Y','RER','T'};
timeline=1947:0.25:2022.25;
if length(timeline)~=length(data_US)
    error("timeline does not match")
else
    data_US=data_US(timeline>=1983 & timeline<2020,:);
    timeline=timeline(timeline>=1983 & timeline<2020);
end
G_detrended=detrend(log(data_US(:,1)));
Y_detrended=detrend(log(data_US(:,2)));
T_detrended=detrend(log(data_US(:,4)));

temp=[T_detrended, G_detrended, Y_detrended, log(data_US(:,3))];
setup.var_names ={'T','G','Y','FX (up is apprecation)'};

YY=temp(p+1:end,:);
XXact=[];
for lag=1:p
    XXact=[XXact temp(p+1-lag:end-lag,:)];
end

n = size(setup.var_names,2); % Number of variables
GDPRATIO =[1 0.2 0.2];

% Fiscal elasticities for Blanchard Perotti identification
BP_EETA_SIMPLE      = zeros(n-2,nP);
BP_EETA_SIMPLE(1,1) = 1.7;

%-------------------------
% Generate Minnesota Prior
%-------------------------

% tau    : Overall tightness
% d      : Scaling down the variance for the coefficients of a distant lag
% w      : Number of observations used for obtaining the prior for the
%          covariance matrix of error terms (usually fixed to 1).
% lambda : Tuning parameter for coefficients for constant.
% mu     : Tuning parameter for the covariance between coefficients*/

%--------------------
% Dummy Observations
%--------------------

tau	   =   0.2;
d	   =   2.5;
w	   =   0;
lambda =   0.5;
mu	   =   0.5;

%get "presample" moments
YY0     =   YY(1:T0,:);
ybar    =   mean(YY0)';
sbar    =   std(YY0)';
pre_sample_moments  =   [ybar sbar];

%------------------------------------------
% Generate matrices with dummy observations
%------------------------------------------

hyp = [tau; d; w; lambda; mu];
[YYdum, XXdum, breakss] = varprior_h(n,p,nex_,hyp,pre_sample_moments);

%--------------------
% Build matrix of actual observations
%--------------------
YYact = YY;
nobs    = size(YY,1);  % number of observations
XXact = [XXact ones(nobs,1)];

%-------------------------
% Estimation Preliminaries: get OLS estimates and construct companion form
%-------------------------

J = [eye(n);repmat(zeros(n),p-1,1)]; % Page 12 RWZ; % J: shock linear combination selector matrix
F = zeros(n*p,n*p);    % Matrix for Companion Form
I  = eye(n);
for lag=1:p-1
    F(lag*n+1:(lag+1)*n,(lag-1)*n+1:lag*n) = I;
end

X = [XXdum; XXact];
Y = [YYdum; YYact];
T = size(X, 1);
ndum = size(XXdum, 1);

% Compute OLS estimates
B = (X'*X)\(X'*Y);       % Point estimates
U = Y-X*B;               % Residuals
Sigmau = U'*U/(T-p*n-1); % Covariance matrix of residuals
F(1:n,1:n*p)    = B(1:n*p,:)';


[F21Ord1BP, IRF_scaled_BP, structural_shocks_BP]= vm_mult_simple(F,J,Horizon,BP_EETA_SIMPLE,U,T,n,p,nP,GDPRATIO); 
[n_periods,n_vars,n_shocks]=size(IRF_scaled_BP);

structural_shocks.BP_T=structural_shocks_BP(end-nobs+1:end,1);
structural_shocks.BP_G=structural_shocks_BP(end-nobs+1:end,2);

[~, A0PF, Ltilde] = penalty_fun_simple(Sigmau,n,p,F,J,Horizon,nperPF);

PF_EETA_SIMPLE      = zeros(n-2,nP);
pos.T=1;
pos.G=2;
pos.Y=3;
shock.BC=1;
shock.G=2;
shock.T=3;

% get elasticity of T and G w.r.t. other variables (Y here)
PF_EETA_SIMPLE(pos.Y-nP,pos.T) = -A0PF(pos.Y,shock.T)/A0PF(pos.T,shock.T); % First row tax, third row GDP; Second column is tax shock;
PF_EETA_SIMPLE(pos.Y-nP,pos.G) = -A0PF(pos.Y,shock.G)/A0PF(pos.G,shock.G); % Second row spending, third row GDP; Third column is spending shock;
[F21Ord1PF, IRF_scaled_PF,structural_shocks_PF] = vm_mult_simple(F,J,Horizon,PF_EETA_SIMPLE,U,T,n,p,nP,GDPRATIO);

[n_periods,n_vars,n_shocks]=size(IRF_scaled_PF);
structural_shocks.PF_T=structural_shocks_PF(end-nobs+1:end,1);
structural_shocks.PF_G=structural_shocks_PF(end-nobs+1:end,2);
timeline=timeline(end-nobs+1:end);
save('results/structural_shocks_point_1983','structural_shocks','YY','timeline');
clear('structural_shocks','A0PF','PF_EETA_SIMPLE','temp','data_mat','Ltilde')

N0=zeros(size(X',1),size(X,2));
nnu0=0;
nnuT = T +nnu0;
NT = N0 + X'*X;
Bbar0=B;
S0=Sigmau;
BbarT = NT\(N0*Bbar0 + (X'*X)*B);
ST = (nnu0/nnuT)*S0 + (T/nnuT)*Sigmau + (1/nnuT)*((B-Bbar0)')*N0*(NT\eye(n*p+nex_))*(X'*X)*(B-Bbar0); %% Constant (check)
STinv = ST\eye(n);

%------------------------------------------------------------
% Define objects that store the draws
%------------------------------------------------------------
BP_DM_T_SIMPLE = zeros(nd-bburn,Horizon);
BP_DM_G_SIMPLE = zeros(nd-bburn,Horizon);
IRF.BP_simple=NaN(Horizon,n,3,nd-bburn);
IRF.BP_simple_balanced_budget=NaN(Horizon,n,nd-bburn);
BP_EETA_TG = zeros(nd-bburn,1);
PF_DM_T_SIMPLE = zeros(nd-bburn,Horizon);
PF_DM_G_SIMPLE = zeros(nd-bburn,Horizon);
PF_EETA_SIMPLE_store = zeros(2,nd-bburn);
IRF.PF_simple=NaN(Horizon,n,3,nd-bburn);

PF_EETA_SIMPLE      = zeros(n-2,nP);
record=0;
counter = 0;

disp('                                                                  ');
disp('        BAYESIAN ESTIMATION OF VAR VIA BLOCK MCMC                 ');
disp('                                                                  ');

%-----------------
% MCMC Chain
%----------------

while record<nd
    %------------------------------------------------------------
    % Gibbs Sampling Algorithm
    %------------------------------------------------------------
    % STEP ONE: Draw from the B, SigmaB | Y
    %------------------------------------------------------------
    % Step 1: Draw from the marginal posterior for Sigmau p(Sigmau|Y,X)
    get_draw =1;
    iter=1;
    R=mvnrnd(zeros(n,1),STinv/nnuT,nnuT)';
    Sigmadraw=(R*R')\eye(n);

    % Step 2: Taking newSigma as given draw for B using a multivariate normal
    bbeta = B(:);
    SigmaB = kron(Sigmadraw,NT\eye(n*p+nex_));
    SigmaB = (SigmaB+SigmaB')/2; %assure symmetry
    while any(abs(eig(F))>1) || get_draw      
        Bdraw = mvnrnd(bbeta,SigmaB);
        % Storing unrestricted draws
    
        Bdraw = reshape(Bdraw,n*p+nex_,n);% Reshape Bdraw from vector to matrix
        Udraw = Y-X*Bdraw;      % Store residuals for IV regressions
    
        Sigmadraw=cov(Udraw);
        Udraw=detrend(Udraw,0);

        F(1:n,1:n*p)    = Bdraw(1:n*p,:)';
        get_draw=0;
        iter=iter+1;
        if iter==1000
           error('Unable to find stable draw')
        end
    end
    
    record=record+1;
    counter = counter +1;
    if counter==0.05*nd
        disp(['         DRAW NUMBER:   ', num2str(record)]);
        disp('                                                                  ');
        disp(['     REMAINING DRAWS:   ', num2str(nd-record)]);
        disp('                                                                  ');
        counter = 0;
    end

    if record > bburn

        if do_BP == 1
            [F21Ord1BP, IRF_scaled, structural_shocks]= vm_mult_simple(F,J,Horizon,BP_EETA_SIMPLE,Udraw,T,n,p,nP,GDPRATIO); 
            BP_DM_T_SIMPLE(record-bburn,:) = F21Ord1BP(:,1);
            BP_DM_G_SIMPLE(record-bburn,:) = F21Ord1BP(:,2);
            IRF.BP_simple(:,:,:,record-bburn)=IRF_scaled;

        end

        if do_PF == 1
            [~, A0PF, ~] = penalty_fun_simple(Sigmadraw,n,p,F,J,Horizon,nperPF);
            PF_EETA_SIMPLE(pos.Y-nP,pos.T) = -A0PF(pos.Y,shock.T)/A0PF(pos.T,shock.T); % First row tax, third row GDP; Second column is tax shock;
            PF_EETA_SIMPLE(pos.Y-nP,pos.G) = -A0PF(pos.Y,shock.G)/A0PF(pos.G,shock.G); % Second row spending, third row GDP; Third column is spending shock;
            PF_EETA_SIMPLE_store(:,record-bburn)=PF_EETA_SIMPLE(1,:);
            [F21Ord1PF, IRF_scaled] = vm_mult_simple(F,J,Horizon,PF_EETA_SIMPLE,Udraw,T,n,p,nP,GDPRATIO);
            PF_DM_T_SIMPLE(record-bburn,:) = F21Ord1PF(:,1);
            PF_DM_G_SIMPLE(record-bburn,:) = F21Ord1PF(:,2);
            IRF.PF_simple(:,:,:,record-bburn)=IRF_scaled;

        end

    end
end

SVAR.BP_DM_T_SIMPLE_Q     = quantile(BP_DM_T_SIMPLE,ptileVEC);
SVAR.BP_DM_G_SIMPLE_Q     = quantile(BP_DM_G_SIMPLE,ptileVEC);
SVAR.BP_EETA_TG_Q         = quantile(BP_EETA_TG,ptileVEC);
SVAR.PF_DM_T_SIMPLE_Q     = quantile(PF_DM_T_SIMPLE,ptileVEC);
SVAR.PF_DM_G_SIMPLE_Q     = quantile(PF_DM_G_SIMPLE,ptileVEC);

PF_EETA_SIMPLE_quantiles=quantile(PF_EETA_SIMPLE_store,ptileVEC,2);

var_names=setup.var_names;
plot_IRFs;

save('results/IRFS_symmetric_1983','BP_G_IRFS',"PF_G_IRFS","var_names","YYact",'ptileVEC')